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ABSTRACT 

We present two-dimensional slab-jet simulations of jets in inhomogeneous media con- 
sisting of a tenuous hot medium populated with a small filling factor by warm, dense 
clouds. The simulations are relevant to the structure and dynamics of sources such as 
Gigahertz Peak Spectrum and Compact Steep Spectrum radio galaxies, High Redshift 
Radio Galaxies and radio galaxies in cooling flows. The jets are disrupted to a degree 
depending upon the filling factor of the clouds. With a small filling factor, the jet re- 
tains some forward momentum but also forms a halo or bubble around the source. At 
larger filling factors channels are formed in the cloud distribution through which the 
jet plasma flows and a hierarchical structure consisting of nested lobes and an outer 
enclosing bubble results. We suggest that the CSS quasar 3C48 is an example of a low 
filling factor jet - interstellar medium interaction whilst M87 may be an example of 
the higher filling factor type of interaction. Jet disruption occurs primarily as a result 
of Kelvin-Helmholtz instabilities driven by turbulence in the radio cocoon not through 
direct jet-cloud interactions, although there are some examples of these. In all radio 
galaxies whose morphology may be the result of jet interactions with an inhomoge- 
neous interstellar medium we expect that the dense clouds will be optically observable 
as a result of radiative shocks driven by the pressure of the radio cocoon. We also 
expect that the radio galaxies will possess faint haloes of radio emitting material well 
beyond the observable jet structure. 

Key words: hydrodynamics - ISM: clouds - galaxies: active - galaxies: ISM - galax- 
ies: jets - galaxies: individual (3C48) 



1 INTRODUCTION 

Previous models of jets propagating through the interstel- 
lar medium of radio galaxies have almost always assumed 
a uniform, usually hot (T ~ 10 7 K) interstellar medium. 
These models have contributed substantially to our under- 
standing of the propagation of jets in evolved radio galaxies 
wherein much of the dense matter that could have previ- 
ously obstructed the jets has been cleared away. However, 
it is now clear that in many young radio galaxies, i.e. Giga- 
hertz Peak Spectrum (GPS) and Compact Steep Spectrum 
(CSS) sources, substantial interaction between the expand- 
ing radio plasma and dense interstellar clouds occurs. This 
interaction has a number of effects: Two of the most impor- 
tant are that it significantly distorts the radio morphology 
and it also gives rise to shock excited line emission. Strong 
observational evidence for the latter effect has been obtained 
from combined Hubble Space Telescope (HST) and radio ob- 
servations of a number of CSS sources ide Vries et al-lfl 999t 



lLabiano et ail 120031: lO'Dea et all l2003h showing a strong 
alignment effect between optical and radio emission and 
clear evidence for disturbance of the optically emitting gas 
by the expanding radio lobe. On the other hand, the dense 
gas cannot have a large filling factor in the case of GPS 
and CSS sources since the expansion spe ed of the radio 
plasma is usually quite high: ~ 0.1 — . 4 c iConwavll2002b ; 
i Murgia , et aT)ll999l: iParma et alJll999t iMurgia et al] I200S : 
iMurgial l2003h . Thus the picture that is conveyed by the 
combined radio and optical observations is of radio plasma 
propagating through a predominantly hot and tenuous in- 
terstellar medium in which dense clouds are embedded with 
a low filling factor. Nevertheless, as we show in this paper, 
even a low filling factor of dense gas can have a pronounced 
effect on the evolution of the radio lobes. 

High redshift radio galaxies provide another example in 
which radio plasma - cloud interactions are clearly impor- 
tant. For example, in the high redshift radio galaxy 4C41.17, 
interactions between the jets and large clouds in the Lyman- 
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a emitting halo of the parent galaxy lead to shock-excited 
line emission and to accelerated rates of star formation 
( Bick nell et, alJ Il998h . Returning to the low-redshift Uni- 
verse, radio plasma - clumpy ISM interactions may be re- 
sponsible for the heating of the ISM that is necessary to 
counteract the otherwise large mass inflow rates that are 
inferred from X-ray observations of cooling flows. 

Thus, there is ample physical motivation for detailed 
investigations of the manner in which jets and the related 
radio lobes interact with realistic interstellar media. One of 
the most significant drivers of this research is that all of the 
scenarios discussed above represent exam ples of AGN-galaxy 
feedback which, it ha s been argued (e.g. ISilk fe Reeslll998t 
iKawakatu et al.ll2003t) . is a significant determinant of the 
ultimate structure of galaxies and in particular, the relation 
between the mass o f the central black ho l e and the mass 
of th e galaxy bulge iM agorrian et al.M 1998: Tre maine et alJ 
120021) . 

In this first paper in this series we explore the nature of 
the interaction between jets and their lobes through a series 
of two-dimensional slab jet simulations. The advantages of 
this approach is that two dimensional simulations offer the 
prospect of high resolution and an initial rapid investigation 
of parameter space, which capture some of the features of 
a more realistic three dimensional simulations. Fully three 
dimensional simulations will be presented in future papers. 

The choice of coordinate system is non-trivial when rep- 
resenting supersonic jets in two spatial dimensions. Many as- 
pects of a jet burrowing through a uniform external medium 
have been adequately described in the past by axisymmetric 
simulations. However an axisymmetric simulation is less ap- 
propriate for physical systems where the jet is expected to 
be bent or deflected. When colliding directly with a dense 
cloud, an axisymmetric jet must either burrow through the 
cloud or cease its advance. However, a two-dimensional slab- 
jet can be deflected, although both jet and cloud are essen- 
tially infinite in the direction perpendicular to the compu- 
tational grid. Our present work is concerned with systems 
where jet-cloud collisions and indirect interactions domi- 
nate. We therefore opt for a Cartesian grid and simulate 
transient slab jets. 



2 CHOICE OF JET AND ISM PARAMETERS 
2.1 Jet energy flux, cloud and ISM density 

The choice of jet par ameters is based on a self-similar model 
for lobe e xpansion feicknell et all il997f) . based upon the 
model by iBegelmanTi^gfiFo . In this model, the lobe is 
driven by a relativistic jet which expands in such a way 
that £ the ratio of hot spot to lobe pressure is constant. A 
value C ~ 2 was inferred by Begelman from simulations by 
iLind et alJ l|l989]) . Let po be the atmospheric mass density 
at a fiducial distance Ro from the core, let 5 be the index 
of the power variation of atmospheric density with radius 
(p = po(R/Ro) ~ S ) and let the jet energy flux be Fe- The 
temperature of the external medium is unimportant in this 
model since it is the inertia of the background medium that 
governs the expansion of the lobe. The advance speed, Db of 
the head of the bowshock as a function of distance R from 



the core is given by: 

v* = 3 i ( F B V /3 (R) (S ~ 2)/3 

c [18(8-<5)7r]£ \P0R0J Vi?o/ 

where the jet energy flux is 10 Fe,4b er g s_1 , n(100 pc) is 
the number density at 100 pc and the numerical values are 
for £ = 2 and 5 = 1.5. Typical expansion speeds of CSS 
and GPS sources are of order (0.05 — 0.4) c with th e largest 
velocities measured for the most compact sources llConwavl 
l2002bl:lMurgia et alll999h . lConwavl (|2002ari has reviewed es- 
timates of number densities ~ 1cm -3 for three GPS sources, 
and so a fiducial value of n(100 pc) = 1 cm -3 is reasonable. 
Compatibility of the expansion speeds estimated from equa- 
tion Q and the observations indicates jet energy fluxes in 
the range of io 46-47 ergs -1 . Jets of this power are also con- 
sistent with a conventional ratio ~ 10~ 12 — 10~ n Hz -1 of 
1.4 GHz monochromatic radio power to jet energy flux and 
a median power at 5 GHz ~ 10 27 - 5 WHz" 1 . Nevertheless, 
in other sources which morphologically resemble CSS radio 
galaxies, such as the less radio powerfu l Luminous Infrared 
Galaxies studied bv lDrake et all <l2003l) . it is apparent that 
the jet powers are somewhat less. Thus, in our simulations, 
we have adopted jet powers of 10 45 and 10 46 e rg s" 1 . 

Hubble Space Telescope observations llO'Dea et alJ 
l2002Ude Vries et alll997Lll99flh show abundant evidence for 
dense clouds of emission-line gas associated with CSS radio 
galaxies and quasars. For a temperature 10 4 T Ci 4 K of the 
dense clouds, a number density n x ~ 1 cm -3 and tempera- 
ture 10 7 T Xi 7 K for the hot interstellar medium, the number 
density in the clouds is n c » 10 3 n x T x , 7 /T c , 4 cm 3 , assuming 
that the initial, undisturbed clouds are in pressure equilib- 
rium with their surroundings. If we adopt such a value for 
the density of the interstellar medium then equation Q im- 
plies that the velocity of advance would be ~ 0.02 c — much 
lower than what is observed. Hence, as discussed in § Q the 
dense clouds must have a small filling factor. Accordingly, 
in the simulations described below, we create clouds of ap- 
proximately the inferred density and investigate the effect 
of various filling factors. 

With the inferred radio source and cloud geometry, the 
line emission orginates from dense, relatively cool clouds 
that are shocked as they are overrun by the e xpanding radio 
source . A model of this typ e was proposed b vlde Vries et all 
Jl999tl and elaborated by iBicknell etaH (I20031) . In both 
treatments, the excess pressure of the advancing radio lobe 
drives shocks into the dense clouds. In the simulations de- 
scribed below, this is seen to occur but some direct jet-cloud 
collisions are also apparent. For radio lobe driven shocks (the 
main source of shocks), the shock velocity is: 




where n c and n x are the cloud and ISM number densities 
respectively. Shocks in the velocity range of 300 — 950 kms - 
are produced for bow-shock advance speeds ~ 0.1c and den- 
sity ratios n c /n x ~ 10 4 — 10 3 . This is in accord with the ob- 
served velocity dispersions and velocity offsets inferred from 
HST spectra of a sample of CSS sources dO'Dea et all2002h . 
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Radiative shocks do not develop instantaneously; the 
time required for a shock to become fully radiative is effec- 
tively the time, ti, required for post-shock gas to cool to 
~ 10 4 K. In order that there be a display of optical emis- 
sion lines, ti must be less than a fraction / ~ 1 of the dy- 
namical time associated with the expansion of the lobe. In 
estimating t& we utilize a cooling function, which estimates 
the non-equilibrium cooling following a high velocity shock, 
and which is derived fro m the results of the one-dime nsional 
MAPPINGS III code Sutherland fc Dopital Ejj) , This 
cooling finction is also used in our simulations described 
below. The condition on ti becomes a condition on the bow- 
shock velocity 



(Ht) < 0.07/ 2 ( 



10 4 cm- 



n c /n x 
10 4 



f-V 



-(3) 



jBicknell et alJl2003t) . This upper limit on Vb is insensitive 
to the various parameters, in particular, /. If condition 
is not satisfied near the head of the radio lobe it may be 
satisfied at the sides where the bow-shock expansion speed 
is about a factor of 2 lower. In this case one would expect 
to see emission from the sides but not the apex of the radio 
source. This condition gives a limiting bow-shock velocity 
Vb ~ 0.1c. Therefore, it would be interesting to see whether 
the more rapidly expanding sources have a lower ratio of 
emission line flux to radio power, and whether the emission 
line flux is distributed in a different manner to that in lower 
velocity sources. These conditions may also be related to 
the causes of the tr ailing optical emission in CSS sources 
ide Vries et alJll999h . There are two possibilities: (a) The 
apex of the bow-shock is non-radiative (if Vb exceeds the 
above limit) or (b) The number of dense clouds per unit 
volume decreases at a few kpc from the core, i.e. beyond the 
normal extent of the narrow line region. We are inclined to 
favour option (b) but (a) should also be kept in mind. 

The simulations capture these characteristics that de- 
pend on the time of the effective onset of radiative shocks in 
the clouds. We see a range of shock conditions varying from 
almost pure adiabatic shocks to strongly radiative. When 
the cooling time of gas becomes shorter than the Courant 
time for the simulation, the integration of the energy equa- 
tion ensures that the radiative shock occupies approximately 
one pixel in which the gas temperature is equal to the prede- 
fined equilibrium temperature (10 4 K in these simulations). 



2.2 Overpressured jets 
density ratio 



Mach number and 



In order to carry out realistic 2D simulations, we wish to 
establish jets with powers per unit transverse length that 
are relevant to the types of radio galaxies mentioned in the 
introduction, in combination with other reasonable values 
of parameters such as jet width and density. If the number 
density of the interstellar medium is ~ 1 cm -3 on 100 pc 
scales, then the jet must be overpressured. We first show 
this in the context of a sub-relativistic jet model and then 
support the proposition with a relativistic jet model. 

Let a nonrelativistic jet have a density r\ times that of 
the external medium, an internal Mach number M, a specific 
heat ratio 7, a pressure £ times the external ISM pressure 
and a cross-sectional area Aj. Further, let p x and v K be the 



density and isothermal sound speed of the ISM. Then the 
jet mass, momentum and energy fluxes are, collecting the 
dimensional quantities on the left: 



M = 

n = 

F E = 



(p x « x Aj) My/lfr 

(pxvlAi) £(1 + 7 M 2 ) 
\ 7- 1 



-y3£3 

V 



(4) 
(5) 

(6) 



Equation @ clearly shows that the energy flux can 
be made arbitrarily large by decreasing the density ratio 
r\. However, there is a natural lower limit on)) - the value at 
which an intially relativistic jet becomes transonic. The crit- 
ical v alue is defined by the condition pj c 2 /4pj = 1 teicknell 
11994ft . implying that Vclit = 4£(« x /c) 2 w 5.9 x 10" 6 ^T x , r . 
With the overpressure factor, f = 1, r] = 10 , a jet diam- 
eter of 10 pc, T x = 10 7 K and M = 20, the energy flux is 
only ~ 10 44 erg s _1 . One requires an overpressure factor of 
about 100 to achieve 10 46 ergs -1 . 

This result is easily verified for a relativistic jet whose 
energy flux is given by 



F E = {p^vlcAi) 4£/3 jet r 2 et 



1 + 



-jet 



1 



(7) 



where \ ~ PjC 2 /4pj, and pj and pj are the jet rest-mass 
density and pressure respectively. If £ = 1, T = 10, a jet 
diameter of 10 pc, T x = 10 7 K and x = (the latter be- 
ing appropriate for a jet dominated by relativistic plasma), 
then the energy flux is only 1.2 x 10 43 erg s _1 . Hence, for a 



10 z 



in order to achieve the 



relativistic jet, we require £ - 
required jet power. 

Is an overpressured jet physical? The overpressure ratio 
refers to the pressure with respect to the undisturbed inter- 
stellar medium. However, the lobe itself is over-pressured as 
a result of the advancing bow shock so that it is possible 
for the jet to be in approximate equilibrium with the ambi- 
ent lobe. In a real galaxy the over pressured lobe may not 
extend back to the core, in which case an overpressured jet 
would behave initially like a classical overpressured, under- 
expanded jet. 

Therefore, in our simulations, we have adopted low val- 
ues of the density ratio, rj ~ 10 -3 , Mach numbers, M of the 
order of 10 and a high overpressure ratio in order to obtain 
the necessarily high jet energy fluxes (or their equivalent 
in 2D). Such a low value of r\ has the additional effect of 
producing superluminal velocities (/3 ~ 5 — 6) in some simu- 
lations, highlighting the shortcomings of simulating power- 
ful jets with a non-relativistic code. However, the essential 
point is that it is mainly the jet power that is physically 
important in jet-cloud interactions. 

One could adopt the strategy of adopting non- 
relativistic parameters that correspond more closely to those 
of relativistic jets, for example by equating the velocity 
an d pressure and matching ei ther the power or thrus t (e.g. 
see lKomissarov fc Falld Jl99rJ) and lRosen et"aH \ 1999ft 1. We 
have not done this explicitly here since there is no trans- 
formation that perfectly maps a non-relativistic flow into a 
relativistic equivalent. Clearly, this initial investigation of 
parameter space needs to be followed by more comprehen- 
sive simulations utilising a fully relativistic code. 



4 C. J. Saxton et al. 



2.3 Prescription for cloud generation 

In our simulations we prescribe a uniform density for the 
hot interstellar medium. Within this medium we establish 
an initial field of dense clouds that mimics the distribution 
of clouds in the interstellar medium, with a random distri- 
bution of positions, shapes and internal density structures. 
We consider that this is more desirable than specifying, say 
circular clouds, albeit with random centres; in that case the 
subsequent evolution may be governed by the initial sym- 
metry. An advantage of the approach that we describe here 
is that the distribution of the sizes, shapes and directions of 
propagation of radiative shocks is more realistic than that 
obtained from initially highly symmetric clouds. 

We prescribe the density and distribution of clouds us- 
ing a Fourier approach. In the first step we establish a grid 
of complex numbers in Fourier space and assign a random 
phase to each point. The largest wave number is compatible 
with the spatial resolution of the simulations. An amplitude 
filter is applied, oc |k|~ 5/3 , in order to mimic the power spec- 
trum of structures formed by turbulence. Taking the Fourier 
transform of the wavenumber data produces an array of den- 
sity fluctuations in real space. These fluctuations, cr(r), are 
normalised so that their mean is and their variance is 1. 
The final number density of gas is then defined by the rela- 
tion: 



iimin + n CT ((j — (J m in) wherever a > a„ 
n x elsewhere 



(8) 



where n x is the number density of hot inter-cloud gas, n m i n 
is the minimum number density of gas in a dense cloud, n CT 
is a scaling factor and a m i n is a threshold value determining 
whether a particular cell is part of a cloud or else hot inter- 
cloud medium. In practice we choose the mean density of 
cloud gas, n and set the value of n a implicitly. 

We apply upper and lower limits to the density distribu- 
tion in k-space. Amplitudes are zeroed for k < fc m i n , where 
fcmin corresponds to the Jeans wave number. This prevents 
the insertion of clouds that would have collapsed to form 
stars and also avoids the generation of a cloud field domi- 
nated by a single large cloud-mass. Amplitudes are also ze- 
roed for k > fc ma x where fc max corresponds to 3 or 4 pixels of 
spatial resolution. Fine features, at the 3-4 pixel level, can 
introduce excessive and artificial sensitivity to effects on the 
level of individual pixels. 

Finally, radial limits are applied to the cloud field as 
a whole. An outer limit, r max , (typically ~ 0.5 — 5 kpc) 
limits the extent of the cloud field. The inner radius, r m i n , 
is defined by the Roche limit in the vicinity of a black hole: 



( 3 M hh 
\4tt p c 



1/3 



(9) 



If the black hole mass Mbh = 10 9 Afo and the cloud density 
p c ~ 15M Q pc~ 3 then r m i n ~ 250 pc. 

Table presents a summary of the properties of the 
clouds and intercloud medium, and the extent of cloud cov- 
erage in our simulations (which are each described in specific 
detail in fJ|J. Figures and [5] represent the initial distribu- 
tions of cloud and intercloud gas in three of our simulations. 
Figures |3] and 0] show how the masses of individual clouds 
are distributed in the initial conditions of four of our simu- 
lations. As the cumulative distribution functions show, most 



clouds have masses of less than a few thousand solar masses. 
The scatter plots show the characteristic depth of clouds: in 
mean terms the distance from an interior point to the sur- 
face, r = DI /S where D is the number of spatial dimensions, 
I is the interior (area of a 2D cloud) and S is the surface 
(perimeter of a 2D cloud). A small r means that most of 
the cloud mass is well exposed to surface effects, such as 
ablation and the penetration of shocks from outside. If r is 
large than the interior regions are typically deeper from the 
locally nearest surface, and correspondingly better insulated 
from effects at the surface. 

Large clouds, near the Jeans limit (e.g. thickness 2r ~ 
105 pc in the B series) are rare in every simulation. The sim- 
ulations with a greater cut-off a m in, and hence greater fill- 
ing factor of clouds, also feature individual clouds of greater 
mass than the simulations with lower o" m i n . 



2.4 Code 

In the simulations reported here, we have used a lo- 
cally modified version of a Piecewise Parabolic Method 
JColella fc Woodwardll984T) code, known as VH-1, published 
on the web by the computational astrophysics groups at the 
Universities of Virginia and South Carolina 1 . Our modifi- 
cations of the VH-1 code include restructuring of the code 
for efficient computation on the computers at the Australian 
National University (formerly a Fujitsu VP300; now a Com- 
paq AlphaServer SC), addition of code for dynamically im- 
portant cooling and other code for coping with a numeri- 
cal shock instability ^Sutherland et all l2003'l. The code has 
already been used in several adiabat ic applications (e.g. 
ISaxton et al jEobH a. ISaxton et al.ll2002l b) and most recently 
in two-dimensional simu lations of fully radiative wall shocks 
Sutherland et al.l Eoolj) . In the latter paper, cooling was 
dynamically important and the code for preventing numer- 
ically unstable shocks was implemented and proved to be 
important in elminating physically unrealistic features. We 
refer to this augmented code as ppmlr. 

We have recently ported the ppmlr code to the Message 
Passing Interface (MPI) environment (see Appendix lAt . 
This port enables us to take advantage of the Compaq Al- 
phasever parallel processor architecture. The MPI port is de- 
scribed in the appendix. The simulations described here typ- 
ically involve 16 processors. Future three-dimensional simu- 
lations will regularly use 256 processors and in exceptional 
circumstances, up to 512 processors. 



3 SIMULATIONS 

In the following subsections we describe a number of simula- 
tions that exhibit a wide range of phenomena when we vary 
fundamental parameters such as the jet energy flux (Fe), 
the jet density ratio (rj) and the filling factor of clouds (/ c ) 
obstructing the path of the jet. The jet parameters in the 
simulations are presented in Table [5] and the characteris- 
tics of the clouds and intercloud media are given in Table 



1 See http:// wonka.physics.ncsu.edu/pub/VH-l/ See also 
Blondin & Lukfin ( 1993) for the extension of the PPM method 
to curvilinear coordinates. 
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Figure 1. Logarithmically scaled map of density relative to the external medium, p/p x , in the initial cloud field of the simulation Al. The 
dark gray region represents the intercloud medium and the lighter gray regions represent cloud cells. Each cloud has an internal density 
structure defined by the Fourier construction described in the text ( 42. 'dt . The region shown extends over 5 kpc X 5 kpc (1600 X 1600 cells), 
with the jet origin at the middle of the left edge. 



Table 1. Parameters and properties of the inter-cloud medium and initial cloud fields in the respective 
simulations. The X-ray emitting inter-cloud medium has temperature T x and particle number density n x . 
Clouds have a minimum density n c m ; n , mean density n c . The cloud generation threshold is <r m ; n and the 
clouds occupy a total area A c , corresponding to a filling factor of / c . The total exposed cloud perimeter is 

Pc 



Model 


Tx 
(K) 


n x 
(cm~ 3 ) 


^c,min/ ™x 


n c /n x 




A c 
(pc 2 ) 


fc 




Pc 
(PC) 


Al 


10 V 


10~ 2 


300 


500 


1.9 


7.28 x 10 b 


2.92 x 10" 


-<i 


5.83 x 10 4 


Bl 


10 7 


1 


500 


1000 


2.6 


1.55 x 10 3 


4.11 x 10" 


-3 


6.16 x 10 2 


B2 


10 7 


1 


500 


1000 


1.9 


1.15 x 10 4 


3.04 x 10" 


-2 


3.51 x 10 3 


B3 


10 7 


1 


500 


1000 


1.2 


4.43 x 10 4 


1.18 x 10" 


-1 


9.74 x 10 4 


B3a 


10 7 


1 


500 


1000 


1.2 


4.43 x 10 4 


1.18 x 10- 


-1 


9.74 x 10 4 
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Figure 3. Initial characteristics of the cloud populations in model Al (c m i n = 1.9). The left column shows a cumulative distribution 
function of the individual cloud masses per parsec. The right column presents a scatter plot of the characteristic depth versus cloud mass 
per parsec. 




10 1 10 2 10 3 10 4 10 5 10 1 10 2 10 3 10 4 10 5 

m (M / pc) m (M / pc) 



Figure 4. Cumulative distribution function of cloud masses (left panel) and mass vs thickness relation (right panel) as in Figure 1^1 but 
for the models Bl (o" m i n = 2.6), B2 (<r m ; n = 1.9) and B3 (cr m i n = 1.2), from top to bottom respectively. 
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Table 2. Parameters and properties of the jets in the respective simulations. M = Mach number, rj = ratio of jet density to intercloud 
medium density, £ is the ratio of the initial jet pressure to the pressure of the external medium, rj is one half the extent of the slab jet (the 
slab jet "radius"). An equivalent three-dimensional jet with circular cross-section has mass flux M, thrust fl and power Fe- 



Model M n f ,\ 



M fl F E m n F E 



(pc) (M Q /yr) (10 33 dyn) (erg/s) A iPx c x A 3 Px cl A jPx 

~~ Al 1Z7 0.384 947 25 3.5 x 10" 2 5A3 io 35 313 2.10 x 10 4 2.10 x 10 s 

Bl 12.7 1.16 x 10~ 3 100 10 1.0 x 10" 2 9.17 10 46 5.58 2.22 x 10 3 1.31 x 10 s 

B2 12.7 1.16 x 10~ 3 100 10 1.0 x 10" 2 9.17 10 46 5.58 2.22 x 10 3 1.31 x 10 s 

B3 12.7 1.16 x 10" 3 100 10 1.0 x 10~ 2 9.17 10 46 5.58 2.22 x 10 3 1.31 x 10 s 

B3a 12.7 1.16 x 10" 3 100 10 1.0 x 10" 2 9.17 10 46 5.58 2.22 x 10 3 1.31 x 10 s 



In each simulation the X-ray emitting inter-cloud medium 
has temperature T x = 10 7 K, and the corresponding isother- 
mal sound speed is our unit of velocity in the simulations, 
= v x = 365 km s . Our unit for time measurement is 
t = v /l pc = 2.68 x 10 3 yr. 

The jet is defined by a group of cells on the left bound- 
ary where the density, pressure and velocity are kept con- 
stant. Elsewhere on the left boundary a reflection condi- 
tion is applied, representing the effect of bilateral symmetry 
of a galaxy with equally powerful jet and counter-jet. The 
boundary cells on other sides of the grid are subject to an 
open boundary condition: in each cycle of the hydrodynamic 
computation, the boundary cell contents are copied into the 
adjacent ghost cells. This one-way boundary allows hydro- 
dynamic waves to propagate off the grid. 

In a journal paper, one is necessarily limited to show- 
ing illustrative snapshots from simulations. Movies of the 
simulations described in the following section, which show 
the detailed dynamical evolution, may be obtained from 
http : //macnab. anu.edu. au/radiojets/ 



3.1 Model Al: high density, overpressured jet 

This simulation exhibits both expected and unexpected fea- 
tures of jets interacting with a distribution of clouds in an 
active nucleus. The jet is initially unimpeded by clouds and 
generally behaves like a jet in a uniform medium showing the 
usual bow-shock, terminal jet-shock and backflow. Then two 
small clouds are overtaken by the bow-shock and the over- 
pressure in the cocoon drives radiative shocks into them. 
Subsequently the jet backflow is distorted by the reverse 
shock emanating from the upper of the two clouds. The in- 
teraction with the cloud field becomes more dramatic when 
the apex of the bow shock strikes one of the larger clouds 
just below the jet axis. Again this drives a strong radiative 
shock into the cloud and the larger reverse shock deflects 
both the jet and backflow sideways. The bow-shock then in- 
tersects the larger of the two clouds and this produces an ad- 
ditional reverse shock that deflects the jet further sideways. 
The jet continues onwards past these two clouds whlist the 
dense clouds in the wake of the bow-shock (i.e. immersed 
within the jet cocoon) are gradually disrupted by the gas 
surrounding them. A large amount of turbulence is gener- 
ated within the jet cocoon by the reverse shocks from the 
bow-shock impacted clouds and the jet flow is clearly af- 
fected but not disrupted. At late times many of the clouds 
immersed within the cocoon show trails resulting from the 
ablation of gas within the turbulent cocoon. However, the 
details of these features will be uncertain until we are able to 



resolve them fully and subdue the effects of numerical mass 
diffusion. 

Notable features of this interaction are: (I) The jet does 
not directly touch any of the clouds. The radiative shocks 
within the clouds are all produced as a result of the over- 
pressure behind the bow-shock. (II) The importance of the 
reverse shocks within the cocoon in deflecting the jet and 
making the cocoon more turbulent. (Ill) The survival of the 
jet despite these interactions (aided by heaviness of the jet 
in this simulation, 77 = 0.3835). 

The left panel of Figure 0shows the jet interacting with 
clouds during an early stage of the simulation, t = 5to. The 
bow shock and therefore the influence of the jet have only 
propagated about 1.7 kpc forward and 0.5 kpc laterally, and 
most of the clouds remain in their initial condition. The 
upper panels of Figure |S| show details of the subregion im- 
mediately surrounding the jet at this time of the simulation. 
Within the high-pressure region bounded by the bow shock, 
all clouds are shocked. The cloud shocks are the intensely 
cooling regions on the cloud perimenters visible in the upper- 
right panel of Figure |5] The cloud density is sufficiently high 
that the shocks have penetrated only a few pc at this stage 
— a small depth compared to a cloud thickness. The other 
noteworthy cooling component visible in the cooling map 
(Figure IHl is the X -ray emission from the bow shock. This 
is more intense near the apex of the bow shock than at its 
sides. 

The lower panels of Figure|S]show the density and cool- 
ing at a later time, t = Wto, when the jet is more advanced 
and the clouds are more eroded. Some clouds that lay in 
the path of the jet at time 5to (located at (0.85, 0.15) kpc 
and (0.95, 0.10) kpc) are now completely shocked and de- 
stroyed. Other clouds have survived closer to the nucleus 
but with less damage because they are out of the jet's path, 
e.g. those at (0.3, 0.4) kpc and (0.3, 0.25) kpc. The jet sur- 
vives an oblique interaction with one group of eroded clouds 
near (0.9, 0.15) kpc, but its width increases to the right of 
this point. A second cloud collision near (1.7, 0.0) kpc de- 
flects the jet into the negative y direction, and causes its 
decollimation into a turbulent bubble of hot plasma accu- 
mulating at the apex of the bow shock. The head of the jet 
is unsteady and changes direction completely within a few 
to. 

The density distribution at time t = 19fo, near the 
finish of the simulation, appears in the right panel of Fig- 
ure Clouds are visible in all stages of shock compression. 
The two major cloud collisions, around (1.00, 0.05) kpc and 
(1.75, 0.00) kpc, now leave the jet well collimated and rel- 
atively straight. The jet loses collimation at a third colli- 
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Figure 2. Density maps of the initial cloud fields of the simula- 
tions Bl, B2 and B3 (from top to bottom respectively). These cloud 
fields were generated from the same fractal distribution but use 
a different threshold <r m i n to define the cloud/intercloud bound- 
aries. The entire grid covers 4 kpc X 4 kpc (2000 X 2000 cells), but 
the square subregion shown above is 1.2 kpc on a side. The jet 
origin is at the middle of the left edge of each panel. 



sion, around (2. 5, ~0. 2) kpc. The turbulent bubbles around 
the head of the jet have grown, and spurt between the clouds 
in several different directions. As was the case at t = 5£o and 
10to, the clouds situated off to the sides have survived much 
better than clouds that lay close to the direct path of the 
jet. 

In summary, whether a cloud behaves radiatively and 
contracts, or adiabatically and ablates is affected by its size, 
location and density. Small clouds are more effectively adi- 
abatic and less durable than large clouds. This may be 
partly because their exposed surfaces are large compared 
to their internal volume, and partly because they are ge- 
ometrically narrower and on average less dense than large 
clouds. Some of the smallest swell and develop ablative tails 
almost as soon as they enter the cocoon, e.g. those between 
(1.3,"0.7) kpc and (2.1, "0.7) kpc when t = 10t , in Figure© 
A larger cloud near (0.2,~0.5) kpc was within the cocoon 
for most of the simulation and yet it was barely affected, 
even until t = 19to (right panel, Figure • Qualitatively, 
each cloud appears to suffer one of three possible fates. 
Small clouds are destroyed quickly regardless of their po- 
sition. Large clouds near the path of the jet affect the jet 
through reflection shocks (even if the jet doesn't strike di- 
rectly) but they inflate and ablate away in the long term. 
Large bystander clouds situated well off the jet axis endure 
almost indefinitely 

Interaction with clouds alters the propagation and 
shape of the bow shock. The bow shock is severely retarded 
when it strikes a cloud and refracts through the gaps be- 
tween clouds. Thus the bow shock surface is conspicuously 
indented ("dimpled") for as long as it remains amidst the 
cloud field. We would expect this effect to manifest in 3D 
simulations also, although the extra degree of freedom in 3D 
is also likely to bring about some differences. 

The jet generally remains identifiable and collimated 
after deflection in two or at most three cloud interactions. 
These interactions need not be direct: it is sufficient for the 
jet to interact with the reverse shocks that reflect off a cloud 
lying within several times r- 3 off the jet's path. These back- 
washing shocks spread continually in the turbulent cocoon; 
they appear to act like a "cushion" against collision between 
jet and cloud. More than three consecutive cloud interac- 
tions disrupts the jet. 

3.2 Low density, overpressured jets 

3.2.1 Jet stability properties, modelBO: 

Our main series of simulations depicts a more powerful jet 
(Fe = 10 46 erg s _1 ) with a much lower density contrast 
(77 = 1.16 x 10~ 3 ) and a large overpressure (£ = 100). The 
initial internal Mach number of the jet is M = 12.72 at the 
nucleus. The jet radius is rj = 10 pc, the square grid is 4 kpc 
wide and the resolution is 2000 x 2000 cells. The background 
gas again has a temperature T x = 10 7 K but the density is 
raised to n x = 1 cm -3 . These parameters are suggested by 
the observation s of GPS sources (e.g. in the summary by 
IConwavll2002ab . 

One of the aims of this work is to distinguish the effects 
of an inhomogeneous medium from the phenomena that are 
intrinsic to the jet itself. Therefore we have performed a 
control simulation containing only the jet and an initially 
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Figure 5. Logarithmically scaled map of p/p x in the model Al at time t = 5io (left panel) and 19fo (right panel). The intercloud medium 
has number density ri x = 0.01 cm -3 . 



uniform intercloud medium. The main large-scale features 
at the end of this simulation (see Figure |7J are: an outer 
shell of highly overpressured gas processed through the bow 
shock, enclosing a turbulent cocoon where jet plasma mixes 
with external gas. The length of the jet is small compared 
to the diameter of the bubble surrounding it. The jet is 
highly unstable: it stays recognisably straight throughout its 
first 0.1 kpc but at greater distances from the origin it flaps 
transversely with increasing amplitude, losing collimation 
at positions beyond ~ 0.5 kpc from the origin. This high- 
lights a feature of two-dimensional simulations that may 
not be present in three dimensions. The instability in this 
jet is driven by 2D turbulence in the cocoon that results 
from the mixing of gas at the edge, although the instability 
is not as rapid as in the ensuing cases with cloudy media. 
The turbulent driving is probably stronger here because of 
the persistence of vortices associated with the well-known 
inverse cascade in two dimensions. The instability evident 
here is seeded by a very small amount of numerical noise in 

the initial data. 

lHardee fc Norman] lll988l) performed a linear stability 
analysis of the instabilities of two-dimensional slab jets. Two 
qualitative kinds of instabilities occur. An unstable pinch- 
ing mode is responsible for diamond shocks at quasi-periodic 
locations along the length of the jet. Sinusoidal instabilities 
are responsible for transverse oscillations, which disrupt and 
ultimately decollimate the jet. According to the analytical 
and numerical results of lHardee fc Normanl ( |l988l) , the max- 
imally growing fundamental mode has wavelength 
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an e-folding length scale of 
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where a c is the adiabatic sound speed in the cocoon imme- 
diately exterior to the jet surface, the parameter a s j — » 1 for 
M 3> 1, and rj, in this case, is the density contrast between 
the jet and its immediate surroundings within the cocoon. 

In simulation BO, the internal Mach number within the 
jet is reduced to typically M » 6, downstream of the first 
diamond shock (see Figure 0. The cocoon gas immediately 
surrounding the jet has a density typically ~ 4pj and adi- 
abatic sound speed a c ~ 800wo, where vo ~ 365 kms -1 is 
the isothermal sound speed of the undisturbed medium out- 
side the bow shock. Upon substitution into equations 1101 
and 1111 . these values imply exponential growth of sinu- 
soidal instabilities on scales I ~ 100 pc and with wavelengths 
A ~ 200 pc. This is consistent with the growing oscillations 
that are apparent in individual frames of our simulations, 
e.g. in Figure [7] From equation 1121 . the implied period of 
the instability is T* < 0.052 to- The simulated jet thrashes 
back and forth transversely within this typical time. 

The sinusoidal transverse instability revealed in this ho- 
mogeneous ISM simulation are similar in wavelength and 
growth time to instabilities manifest in the inhomogeneous 
simulations discussed below. However, the appearance of 
this instability in this "clean" simulation shows that the 
instability can be properly regarded as a Kelvin-Helmholtz 
instability of the jet itself — albeit, in the inhomogeneous 
cases, stimulated by the turbulence in the cloudy medium. 
Note that growth and wavelength of the transverse KH in- 
stability at first sight seems to be sugge stive of a transonic 
jet. However, the comparison with the lHardee fc Normanl 
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Figure 6. Subrogions of model Al showing the area disturbed by the jet at the times t = 5to (upper row) and lOtrj (lower row). The left 
panels show density p/p x and the right panels show velocity vectors superimposed on maps of the volumetric rate of cooling. The cooling 
is measured in units of erg cm -3 s _1 , but the upper end of the scale is deliberately saturated to keep the bow shock visible. The actual 
maximum rate of radiative cooling in the cloud shocks is 1.8 X 10 3 erg cm~ 3 s" 1 when t = 5tQ. Reflected and refracted shocks are visible 
at distances of order 10 pc around some of the clouds. 



( 1988) expressions for maximally growing wavelength and 
growth time shows that our simultions are consistent with 
the expected behavior of supersonic jets with the quoted 
Mach numbers. 



3.2.2 Small cloud filling factor: model Bl 



This is our simulation with the fewest clouds: their filling 



factor is f c = 4.1 x 10 (o n 



2.6) within a region span- 



ning radially between O.lkpc and 0.5kpc from the nucleus. 
All of the clouds are quickly overtaken by the expanding bow 
shock. As in model Al, each cloud is shocked by the over- 
pressure of the cocoon (up to 5 x 10 4 times greater the ex- 



Interactions of Jets with Inhomogeneous Cloudy Media 11 





ternal pressure), and contracts as it cools radiatively. Some 
clouds are accelerated away from their initial positions by 
the thrust of the jet. For some other clouds, radiative cool- 
ing eventually proves inadequate; the tendency to contract 
is overcome, and the cloud's fate is expansion and abla- 
tion. Halfway through this simulation at time t = 24to, (top 
panel, Figure |Sj, shows the destruction of one such cloud 
in the path of the jet, situated at (0.6, 0.0) kpc but ablat- 
ing away across about 0.2 kpc to the right. Meanwhile, a 
few extremely dense and compact clouds persist in appar- 
ent safety off to the sides of the jet, within a few hundred 
pc of the nucleus. They appear to survive indefinitely and 
appear scarcely changed by the time the simulation ends, as 
seen in the middle panel of Figure[S] In comparison with Fig- 
urcQwe see much more dense gas (derived from the ablated 
clouds) intermixed with jet plasma in the cocoon. 

The last panel in Figure |H| represents the appearance 
of the jet in synchrotron emission at late stages of the 
simulation. Here we h ave used the approximation (e.g. 
ISaxton et alj200ll . l2002Tl that the emissivity j v cx ipp { ' i+a)/2 , 
where p is the pressure, ip is the relative concentration of 
mass originating from the jet (a scalar tracer that is pas- 
sively advected in our code) and a is the spectral index that 
we take to be 0.6. Several diamond shocks are visible in the 
well collimated section of the jet closest to the nucleus. The 
brighter, barred "fishbone" shaped pattern seen at distances 
x > 0.4 kpc is a transient and typical manifestation of the 
transverse instability of the jet (as in § 13. 2. II . 
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Figure 7. Snapshots of the cloudless test model B0 at time 
t = 48to. The upper panels show the local Mach number: the full 
computational grid, and a subregion around the jet. The lower 
panel maps density, p/p x , overlain with velocity vectors. 



3.2.3 Intermediate cloud filling factor: model B2 

In this case the filling factor of clouds is higher with f c = 
3.0 x 10~ 2 , corresponding to cr m m = 1.9. In the first few 
to of simulation time the jet is completely disrupted within 
< 0.1 kpc of the nucleus as a result of the KH instabil- 
ity induced by waves reflected off the innermost clouds. By 
t ~ 12to (see top panel, Figure |UJ, the jet has ablated the 
nearest clouds and cleared an unobstructed region out to 
~ 0.1 kpc on either side. Within this region the jet remains 
well collimated during the rest of the simulation. The clouds 
that define the sides of this sleeve are shocked but they sur- 
vive until the last frame (middle panel, Figure 0. The flow 
of plasma from the head of the jet splits into a fork on either 
side of a cloud group that sits nearly on the jet axis at about 
x ~ 0.5 kpc. Throughout the simulation these clouds re- 
tain dense, efficiently cooling cores but they continually lose 
dense gas entraining into the stream of jet plasma passing 
to either side. Near the end of the simulation (t = 48io, Fig- 
ure |UJ these obstructive clouds are trailed by dense ablated 
gas extending throughout a triangular area up to 0.5 kpc 
further out in the positive x direction. 

The radio emissivity (last panel, Figure |UJ has some 
similarities to the appearance of model Bl, with a bright jet 
remaining well collimated out to around 0.5 kpc. However 
the turbulent plumes of jet-derived plasma beyond the jet's 
head are more disrupted by the dense obstructions of the 
cloud field, Radio bright plasma divides in several directions, 
as Figure|U]shows. The brightest areas of the jet and post-jet 
flows are typically concentrated closer to the nucleus than 
in the less cloudy model Bl. 
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Figure 8. Snapshots of model Bl. Top-left: p/p x and velocity 
field in a subregion at time t = 24i(j. Top-right: t = 48fo, p/px in 
the full grid. Bottom-left: t = 48to, simulated radio emissivity in 
a subregion (the lowest background intensity is 1/ 1024 times the 
peak). 



Figure 9. Snapshots of model B2. Top-left: p/p x and flow ve- 
locities near the jet when t = 12to. Top-right: t = 48to, p/px 
throughout the full grid. Bottom-left: t = 48io, simulated radio 
emissivity around the jet (with a dynamic range of 1024). 
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3.2.4 Large filling factor: model B3 

This is our simulation with largest filling factor of clouds, 
f c — 0.12 and a m in = 1.2. As a result of the high coverage of 
clouds, the jet is disrupted when it is less than 0.1 kpc from 
its origin. Jet plasma streams through the available gaps 
between the persistent clouds (see the density distribution 
in right panel of Figure II It . These low-density, radio-bright 
channels evolve as the cloud distribution changes. A chan- 
nel may close when nearby clouds expand ablatively or are 
blown into positions where they create a new obstruction. 
Turbulent eddies of jet plasma accumulate in pockets within 
the cloud field, most conspicuously in places where a chan- 
nel has closed. The density image, Figure [TTJ1 illustrates the 
initial development of intercloud channels. 

The left panel of Figure fTTI shows simulated radio emis- 
sivity in the inner regions of the galaxy at the end of the sim- 
ulation. The radio-bright material is concentrated near the 
nucleus within an irregular, branched region that is not di- 
rectly related to the alignment of the initial jet. This partic- 
ular configuration of channels is ephemeral, occurring after 
the opening and closure of previous channels as the cloud- 
field evolves. 

Eventually, one or more channels pushes into the open 
region outside the cloudy zone. A spurt of turbulent jet 
plasma accumulates at that point. In later stages of our sim- 
ulations, these eruptions eventually become large compared 
to the cloudy zone, and merge with adjacent eruptions to 
form a diffuse, turbulent, approximately isotropic outer co- 
coon composed of a mixture of jet plasma, ablated cloud gas 
and intercloud gas processed through the bow shock. 

3.2.5 Effect of cooling: model B3a 

We ran a further simulation with the same cloud field as in 
B3 but with radiative cooling switched off (see Figuretl2)l. In 
the absence of radiative cooling, the compression of clouds 
is insignificant. The outflow of jet plasma rapidly turns the 
clouds into ablating, flocculent banana shapes. Large clouds 
and small clouds ablate efficiently. The survival time of indi- 
vidual clouds is shorter than in the simulation including the 
effects of radiative cooling. Dense material is more readily 
ablated from the clouds, and this may block channels of jet 
plasma more readily or at an earlier time, but otherwise the 
radio morphology is qualitatively similar to that of B3. 

3.3 Evolution of the outer part of the radio source 

In each case, as a consequence of the disruption of the jet, 
its momentum is distributed throughout the cocoon in an 
isotropic way. Hence, the long term evolution is similar to 
that of an energy-driven bubble. Figure 1131 presents plots 
of the average radius of the furthermost extent of the ra- 
dio plasma versus time, for represenative simulations. The 
plots also show the corresponding analytic solution for a 
two dimensional bubble with the same parameters of energy 
flux per unit length and density (grey line, see Appendix IBl 
for a derivation). It is evident that the analytic slope of 
the radius/time curves is replicated very well by the simu- 
lations. The offset may be attributed to the time taken for 
the isotropic phase to be established. This may be related to 
the finite minimum width of the jet at its nozzle, the size of 



the innermost nuclear region that is initially cloudless, and 
the intrinsic length scale of the jet's instability (§ 13. 2. II . 

We expect that a similar power-law expansion will occur 
in three dimensions, although this of course is yet to be con- 
firmed. However, given that jet momentum and power can be 
redistributed isotropically in three dimensions, it is reason- 
able to expect that the long term evolution of the outer bub- 
ble will be described by a three-dimensional energy-driven 
bubble. 



4 COMPARISON OF MODEL Bl WITH 3C48 

Whilst these simulations have been two-dimensional, some 
of their features should persist in three-dimensional sim- 
ulations with similar parameters. In particular the trends 
with filling factor will probably persist, to some extent, so 
that comparisons with existing observationsla are of interest, 
even at this stage. We therefore make a tentative compar- 
ison of one of these simulations (model Bl) with the CSS 
quasar 3C48. We can also draw some general conclusions 
from these simulations, concerning the likely overall mor- 
phology of young radio sources. 

Simulation Bl has the lowest filling factor of dense 
clouds and impedes the jet the least. The jet is still disrupted 
but the flow associated with it continues with reduced colli- 
mation and with its forward momentum spread over a much 
larger region. In the disrupted section, transverse oscillations 
are produced, with shocks created at points where the jet 
bends significantly, leading to local enhancements in emis- 
sivity. A large scale bubble is produced outside the cloudy 
region but the jet also retains some coherency in the interme- 
diate region. This particular simulation bears a strong qual- 
itative resemblance to 3C48. In Figure [TH a snapshot of the 
simulated radio emissivity at a particular t ime is compared 
with the 1.6 GHz VLBI image of 3C48 (Wil kinson et alJ 
Il99ll) . In two dimensions, the alternating direction of the jet 
results from the transverse instability that we have already 
discussed ( t!3.2.H . In three dimensions, the corresponding 
instability would be a helical m — 1 mode. In both cases the 
alternating direction of the supersonic flow produces shocks. 
We see evidence for these, in both the snapshot of the emis- 
sivity from the 2D simulation and the image of the jet. 

This particular simulation predicts that there will be 
a diffuse halo of emission external to the jet region of the 
source. This was, in fact, detected in larger scale images o f 
3C 48. One such image presented in I Wilkinson et al] lll99ll) . 
shows diffuse emission on a 1 arcsecond scale in 3C48 ex- 
tending well beyond the compact jet structure evident in 
the VLBI image . More rece nt higher quality images have 
been obtained bv lOihal |2004) and one of their contour maps 
is reproduced in Figure [fTI These images clearly re veal the 
extended structure of 3C48. Wil kinson et alJ <I1991T) had al- 
ready concluded that the morphology of 3C48 is likely the 
result of jet-cloud interactions and cited evidence for asym- 
metric emission line gas northwest of the core as well as the 
presence of dust revealed by the high IRAS flux. They in- 
terpreted the radio knot close to the core as the result of 
a jet-cloud interaction. Whilst this is possible, our simula- 
tion shows that the sort of jet disruption observed in 3C48 
can occur without direct jet-cloud encounters. The general 
turbulence in the jet cocoon can excite jet instabilities, pro- 
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Figure 10. Logarithmically scaled density images (p/p x ) of a central subregion of simulation B3 during the initial erosion of channels of 
jet plasma, and the earliest eruptions of radio bubbles beyond the surface of the cloudy zone. The time steps are t = 3to,6to,9to, 12to in 
the top-left, top-right, bottom-left and bottom-right frames respectively. 



ducing shocks in the supersonic flow either as a result of 
the jet being bent through an angle greater than the Mach 
angle, or though pinching-mode instabilities. (Both of these 
occur in the simulations.) An explanation for the first bright 
knot involving cocoon turbulence may therefore be more fa- 
vorable since direct jet-cloud encount ers are usually highly 
disruptive and I Wilkinson et al.l lll99ll) appealed to a glanc- 
ing deflection in order to explain the radio structure. Given 
that the sort of jet-ISM interaction revealed in our simu- 



lations, is occurring, we expect that the emission line gas 
would show evidence for shock excitation as a result of ra- 
diative shocks i nduce d by the overpressured radio cocoon. 
IChatzichristoul (l200lT) has suggested that the emission line 
spectrum of 3C48 may show evidence for a mixture of AGN 
and shock excitation. 



I Wilkinson et al.l il99lf) also noted out that 3C48, clas- 
sified as a CSS quasar, is much less extended than other 
similar power radio-loud quasars. They speculated that the 
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Figure 11. The final condition of model B3 (t = 48to), comparable to the lower panels of Figures and The left panel is a simulated 
radio emissivity image of the jet's vicinity, with dynamic range 1024 as in previous figures. The right panel represents p/p x throughout the 
full grid. 
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Figure 12. Demonstration of the effects of radiative cooling upon cloud survival. Density frames (p/p x ) are compared at time t = 12to, 
with the left panel showing simulation B3 and the right panel shows simulation B3a, which is the equivalent simulation with radiative 
cooling switched off. In the absence of radiative cooling, clouds expand and ablate rapidly as soon as they pass through the bow shock. 
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Figure 13. Evolution of the bow shock's mean radius compared with the radial growth of an equivalent energy-driven bubble model. The 
left panel shows results for the simulation Al; the right panel shows results for model B3. In each case the grey line is the ideal model, and 
the black line shows simulation results. The dotted line in the right panel shows the cloudless control model BO. 



confinement results from the emission line gas and dust in- 
ferred to be present in the ISM. Qualitatively, we agree with 
this, except to point out that the dispersion of the momen- 
tum flux (by a small number of clouds) prevents a radio 
source from expanding rapidly. As we have shown, the source 
becomes an expanding bubble with a much slower rate of ex- 
pansion than a jet with the same kinetic power. Hence the 
gas and dust is not required to physically confine the radio 
plasma. The slower expansion can be effected by a small 
filling factor of dense clouds that redistribute the jet mo- 
mentum isotropically. 

We also note a qualitative resemblance between the sim- 
ulation B3 and M87. The channeling of the flow in B3 pro- 
duces a hierarchy of structures resembling the hierarchy of 
radio s tructures in M87 apparent in the low frequency radio 
image llOwen et all200oll . However, we have not emphasized 
this with a detailed comparison since the scale of the source 
and the simulation are quite different. Nevertheless, the pro- 
duction of the structure observed in simulation B3 depends 
mainly on the filling factor of the clouds and this mechanism 
should also work on larger scales. We propose to confirm this 
in future with a fully relativistic 3D simulation. 

Other sources with similar morphology to 3C48 and 
which may be being disrupted in a similar way include the 
TeV 7-ray blazars, MKN 501 and MKN 421. The most re- 
cent VLB I images of these sources (Giovannini et alJl2000t 
IPiner et al.lll99sh show similar structures in the form of dis- 
rupted jets and large scale plumes suggesting that these 
sources are disrupted on the sub-parsec scale by a cloudy 
medium. Again the scale of the simulation Bl is such that 
it does not strictly apply. However, the physical principle of 
these jets being disrupted by a small filling factor of dense 
clouds is probably still valid with interestin g repercussions 
for t he production of the 7-ray emission feicknell et alJ 




5 DISCUSSION 

We have presented a series of two-dimensional simulations 
of slab jets propagating through inhomogeneous media. De- 



spite the limitations of the restrictions to two dimensions, 
and to non-relativistic flow, these simulations provide useful 
insight into the characteristics of jets in radio galaxies which 
have an inhomogeneous interstellar medium. These include 
GPS and CSS radio sources, high redshift radio galaxies and 
galaxies at the center of cooling flows. In all of these galax- 
ies there is ample evidence for dense gas clouds that may 
interact with the radio source producing at least some of 
the features evident in our simulations. In the case of GPS 
and CSS sources the inhomogeneous medium may be the 
result of a fairly recent merger; in the case of high redshift 
radio galaxies an inhomogeneous medium is likely to be the 
result of the galaxy formation process at those epochs, as 
large galaxies are assembled from sub-galaxian fragments. 

As well as the limitations to two dimensions and non- 
relativistic flow, there is another feature of our simulations 
that requires some comment and that is the ablation of cloud 
material that is evident in most cases. There is likely to be 
a large contribution of numerical mass diffusion due to the 
very high density contrasts between jet and cloud gas. Al- 
though this has some effect on the evolution, the dominant 
processes in the case of a low filling factor are waves emanat- 
ing from the clouds impacted by the high pressure cocoon 
and, in the high filling factor simulations, direct obstruction 
of the flow plus wave effects. 

Despite the above-mentioned caveats, we have been 
able to deduce these general features from our simulations: 

(i) Jet instability. When propagating through the type 
of inhomogeneous interstellar medium that we have im- 
posed, the jets become unstable. However, this instability 
does not always result from direct jet-cloud collisions but 
occurs indirectly from the waves generated by the impact of 
the high-pressure jet cocoon on the dense clouds. This leads 
to a general level of turbulence in the medium that induces 
unstable Kelvin-Helmholtz oscillations in the jet. An unsta- 
ble jet was also produced in the homogeneous simulation 
(BO) and this was related to turbulence caused by the en- 
trainment of ISM into the cocoon as well as the higher level 
of large scale turbulence in two dimensions. However, the 
simulations with an inhomogeneous medium became unsta- 
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Figure 14. The radio features of 3C48 compared wi th simulatio n Bl. The left panel is a 22 GHz (K-band) VLA A-array image illustrating 
the large-scale (arcsecond scale) structure of 3C48 joiha 2004). This image shows diffuse emission on a larger scale to that of the jet, 
consistent with the prediction of the simulation. The middle panel is a 1.6 GHz VLBI image featuring the jet and bright knots (Wilkinson 
ct al. 1991). At the redshift of 3C48 (z = 0.3670) an angular scale of 0.1 arcsec corresponds to 0.514 kpc (assuming Hq = 70 km s — 1 Mpc -1 , 
fim = 0.27 and f^A = 0.73). The right panel shows a snapshot of radio emissivity from the Bl simulation, at time t = 48<o- Left and middle 
images are reproduced with the kind permission of R. Ojha and P. Wilkinson respectively. 



ble in a shorter distance and we attribute this to the effect 
of the clouds. 

(ii) Production of halo/bubble. The instability of the 
jet is associated with the spreading of the jet's momentum 
and a turbulent diffusion of the jet material into a halo 
or bubble. These are characteristics that we expect to be 
present in three dimensions. The radius-time curve for the 
outer part of the halo is well-described by the analytic ex- 
pression for a two-dimensional bubble. In three dimensions 
a energy driven bubble model should also be appropriate. 
This is useful since it gives us a way of estimating the energy 
flux fed into a large scale radio source from the dimensions 
of the outer halo and the density of the surrounding inter- 
stellar medium. The production of a bubble means that the 
overall increase of size of the radio source with time is dimin- 
ished with respect to a source that propagates unimpeded. 
However, this restriction on the size of the source is not due 
to direct physical confinement of the radio plasma by denser 
gas. Rather, it is the result of the isotropic spread of the jet 
momentum by the process of jet instability described above. 

(iii) The effect of filling factor. The main parameter 
that we have varied in these simulations is the filling factor 
of the dense clouds and we can see that this makes a con- 
siderable difference to the morphology of the radio source. 
With only a small filling factor of clouds, the jet is readily 
disrupted, but retains some directed momentum and pro- 
duces the 3C48-type morphology described in § 21 As the 
filling factor increases, the jet propagates via a number of 
different channels which open and close leading to an hier- 
archical structure of radio plasma that we noted is similar 
to the large-scale structure of M87. 

(iv) The effect of jet density ratio. We also varied the 



jet density ratio parameter, rj, but only in one simulation. 
The simulation Al had r] = 0.38 whereas in the other simula- 
tions rj ~ 10~ 3 . The main noticeable effect is that the heav- 
ier jet is affected less by instabilities and direct jet-cloud 
collisions than the lighter jets. This points to an important 
difference between models with 77 > 0.01 and models with 
rj< 0.001. Of course, such low density ratios are important 
to provide the relevant power for radio galaxy jets. However, 
reassuringly for physics but difficult for modelers, such low 
values of 77 push models into the relativistic regime, as we 
have discussed in § 12.21 There are two important points to 
be made here. The first is that many jet simulations usually 
involve a modest value of r\ > 0.01 and these may be mis- 
leading. Second, the most realistic simulations will be those 
that directly include relativistic physics and which also in- 
corporate the mixing of relativistic and thermal plasma. 

This paper complements other work on the physic s 
of direct jet-clou d collisions (e.g. iFragile et alJ |.2004), 
IWiita et alJ <|200 21) in which the focus has been on the ef- 
fect of a jet on a single cloud. Direct jet-cloud collisions 
are clearly importa nt in the context of so urces such as 
Minkowksi's object (IVan Breueel et allll985l) . In this paper 
we have clearly been more concerned with the overall radio 
source morphology, the evolution of an ensemble of clouds, 
the effects of indirect jet-cloud collisions (i.e. the effect of 
the high pressured cocoon on dense clouds) and trends with 
cloud filling fac tor. However, we note an interesting con- 
nection between IWiita et al] <|2002T) and this work: We have 
both noted that less jet disruption occurs in the case of slow 
dense jets. We also note that the morpholog y of radiative 
cloud shocks deduced by [Fragile et all J2004I) is also repro- 
duced in some of the individual clouds in our simulations. 
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Our simulations have clearly all been two-dimensional. 
What are the expected differences between two-dimensional 
and three-dimensional work? There are at least two (1) In 
two dimensions, vortices survive for a lot longer as a result 
of the inverse cascade of turbulent energy. These vortices 
play a significant role in disrupting the jet and clouds. We 
expect much less disruption from vortical motions in three 
dimensions. (2) The number of independent directions that 
the jet can take in three dimensions is greater so that we 
expect that the clouds will impede the jet less severely and 
that escape from the cloudy region surrounding the nucleus 
will be easier. These issues as well as the inclusion of special 
relativistic effects will be addressed in future papers. 
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APPENDIX A: GRID SUBDIVISION FOR 
MULTIPROCESSOR COMPUTATION 

We confront two practical obstacles to the calculation of 
high resolution hydrodynamics simulations in higher dimen- 
sions: the duration of the run in real-time (which inhibits 
the exploration of a wide parameter space), and the finite 
computer memory (which directly limits the number of grid 
cells that can be computed). We circumvent these limita- 
tions by subdividing the spatial grid into N = 2 n subgrids 
(with n an integer) and assigning each to a program running 
in parallel on a separate processor. The processors communi- 
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cate between themselves using functions of the MPI (Message 
Passing Interface) library 2 . 

If the global grid has dimensions of I x J cells in the x 
and y directions, then we define a "long subgrid" as an area 
with / x (J/N) cells (see the left panel of Figure |B"T1 . and a 
"wide subgrid" as an area of (I/N) x J cells (see the right 
panel of Figure |B"T1 . Within the memory at each processor 
we maintain representations of one long subgrid or one short 
subgrid: see for example the shaded long and wide subgrids 
held by processor 3 in Figure |B"T1 The total number of cells 
in a long subgrid is the same as that of a wide subgrid. 
In fact for each hydrodynamic variable on each processor, 
we identify the long and wide subgrid arrays with the same 
section of memory, which is feasible and reliable since the 
program uses only one of the subgrid representations at any 
time (as we shall detail below). 

We then conceptually subdivide each long or wide sub- 
grid into N blocks separated by transverse cuts (e.g. the 
dotted line divisions in processor 3's subgrids depicted in 
Figure lB"Tll . For each processor p the block numbered p is 
the only area of overlap between the long and wide subgrids 
that it is assigned. 

At the start of the simulation each processor is assigned 
long and wide subgrids according to its identity number. 
These regions are filled according to the initial conditions 
which are defined in term s of global grid coordinates . The 
first hydrodynamic sweep llColella fc Woodw ard 1984.) is in 
the x direction, and is performed by each processor within 
its long subgrid arrays (see left panel of Figure lB"2t . Thus 
the z-sweeps do not encounter any of the subgrid interfaces, 
obviating the need for any communication of flux informa- 
tion between processors at this stage. The one-dimensional 
array vectors associated with each sweeep are the maximum 
possible length, which makes optimal use of processor cache. 

After the x sweep is complete, the evolved density, ve- 
locity, pressure and other zone data are communicated be- 
tween the processors. By use of the MPI_ALLTOALL function, 
the data of block i of the long subgrid of processor j is 
copied to block j of the wide subgrid of processor i, for all 
i and j. These are corresponding areas of the global com- 
putational grid. As the communication sequence proceeds, 
each processor receives data blocks and assembles them into 
complete, ordered arrays describing the wide subgrid that it 
is conceptually assigned. The wide subgrid arrays overwrite 
the memory areas that was previously used for long subgrid 
data. However the sequence of the MPI_ALLTOALL communi- 
cation provides for a complete redistribution of data across 
the set of processors without any collective loss of informa- 
tion. 

When the communication is complete, each processor 
performs a hydrodynamic sweep in the y direction within 
its respective wide subgrid arrays. (See right panel of Fig- 
ure IB2I ') Again, as in the x sweep, no processor bound- 
aries are encountered during the computation, and the one- 
dimensional sweep arrays are long. At the end of the y sweep, 
the processors communicate blocks of their evolved wide grid 
arrays in another MPI_ALLTOALL operation, ensuring that the 
long subgrid arrays are up to date. At this point program 
control operations are performed, e.g. generating output files 
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recording the hydrodynamic variable arrays of each proces- 
sor. 

Care is needed to ensure that the processors evolve 
their subgrids synchronously in simulation time. Within 
each subgrid g the maximum allowable timestep of integra- 
tion for each x or y sweep is determined by a Courant con- 
dition: St g must be less than the characteristic timescales 
of fluid motion or sound waves crossing a cell, e.g. 8t g < 
0.6mm(Sx/v x j,8x/c St i), for all cells i with side length Sx. 
In order to synchronise the simulation across all the sub- 
grids, the processors are made to inter-communicate their 
locally calculated values of St g , and then adopt the global 
minimum allowable timestep, St = min(<5to, Sti, ...Stu-i)- 

This straightforward scheme for multiprocessor PPM hy- 
drodynamics simulations is readily generalised to a three- 
dimensional grid. The 3D analogues of long and wide sub- 
grids are "flat slabs" (full size in the x and y directions, but 
subdivided into N segments in the z direction) and "tall 
slabs" (full size in the x and z directions, but subdivided into 
N parts in the y direction). Using 32, 64, 128 or more of the 
processors of the Australian National University's Compaq 
AlphaServer SC, we have been able to perform 3D simula- 
tions with cubic dimensions of 256, 512 or more cells per 
side (work in preparation). 



APPENDIX B: SELF-SIMILAR SOLUTION FOR 
INFLATION OF A TWO-DIMENSIONAL 
BUBBLE 

Bl Dimensional analysis 

This appendix provides an analytical solution that can be 
used for the comparison of slab-jet simulations with ex- 
pected asymptotic behaviour. The two parameters that de- 
termine the evolution of the bubble are A, the energy flux 
per unit transverse length, and po the density of the back- 
ground medium. In terms of jet parameters, 
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where D is the width of the jet. 

Consider the parameter A/ po. The dimensions of this 
parameter are given by: 
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Let R(t) be the bubble radius at time t. The dimensions of 
Rt~ 3/Ii are also LT ,_3//4 . Hence the self-similar evolution of 
the bubble will be given by: 
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where C is a numerical constant to be determined. The ra- 
dius of the bubble is therefore given as a function of time 
by: 

1/4 
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This derivation is justified and elaborated in the next sec- 
tion, where the constant C is derived. 
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B2 Temporal evolution of bubble 

B2. 1 Formulation of equations 

Consider a cylindrical bubble that is being fed by a jet with 
an energy density per unit length (transverse to the jet) 
equal to A. We consider the jet energy to be effectively ther- 
malised so that the pressure derived from the shocked jet 
plasma drives cylindrical expansion. Let the energy density 
of the bubble be e and the pressure of the bubble be p. 

We represent the bubble in two dimensions in the x — y 
plane and take a length, L of bubble in the z— direction. 
The total energy, E — eirR 2 L of this length of bubble with 
volume itR 2 L is governed by the equation: 



(B5) 



The second term represents the work done by the pressure 
during the expansion; the term on the right represents the 
rate of energy fed into the bubble. Note that we are assumig 
that the dominant form of energy in this bubble is internal 
and that kinetic energy does not contribute significantly. 

Taking out the common factor of L, the evolution of the 
bubble is governed by the equation 

R 2 )+p±(nR 2 ) =A (B6) 
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We solve this equation by introducing the equation of state 
[e = p/(7 — 1)] and by relating the pressure in the bubble to 
its rate of expansion. 

Assuming that the bubble is expanding supersonically 
with respect to the background medium, then its structure 
consists of a strong shock outside the bubble and a con- 
tact discontinuity separating the bubble material from the 
shocked ISM. In this region, the pressure of the shocked ISM 
and the bubble are equal. Suppose that the strong shock is 
moving with velocity v B ^. Then the gas behind the shock is 
moving with speed 3i> s h/4 so that 

R = \ « S h (B7) 

Let the density of the background medium be po- Then the 
pressure behind the shock is given by 
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Note that we have not specified the 7 for the bubble but we 
have assumed that 7 = 5/3 for the background ISM. 

Substituting the equation of state into l|B6fl . and ex- 
panding the first term, the bubble equation becomes: 



7 dR 1 _, 2 dp , 

— —r^P—r- + -tvR 2 -^=A (B9) 

■y-l dtj-ldt y ' 

Then, using equation 1B8I for the pressure, we obtain, after 
some simplification: 

RR 3 + ±R 2 RR= TZI JH. (BIO) 



where a and a are constants. Subsitution into equa- 
tion fMTol gives: 
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Equating the powers of t on both sides: 
3 

a — — 
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as expected. Equating the constants: 
7-1 A 11/4 
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The form of the dependence on the parameter A/po was 
anticipated in ilBll For 7 = 5/3, as used in our simulations: 
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and the solution for R is: 
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B3 Radius of shock wave 

In comparing the analytic solution with simulations it is eas- 
iest to locate the position of the discontunity in the pressure. 
This is the location of the outer shock wave, whose velocity 
is derived from equation (IB7L i.e. 

v sh = I R (B17) 
Hence the radius of the pressure discontinuity is given 



by 
Rp 

for 7 = 5/3 
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This paper has been typeset from a TfrjX/ file prepared 

by the author. 
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B2.2 Solution 

It would be reasonably straightforward to develop a general 
solution of equation 1B10II . However, if we settle for a self- 
similar solution, for which: 
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CPU D CPU 1 CtVl CPU J CPU 4 CPU 5 CPU 6 CPU 7 



Figure Bl. Distribution of the long and the wide subgrid partitions in a simulation run with 8 processors. For illustrative clarification, 
the long subgrid assigned to processor 3 is highlighted in the left panel, and the wide subgrid assigned to processor 3 is highlighted in the 
right panel. Each processor's memory needs only to retain one long subgrid array and one wide subgrid array, plus scalar global variables. 
Collective maps of the global computational grid need only be assembled from output data files after the simulation is run. 
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Figure B2. Hydrodynamic sweeps in the x direction are only performed in the long subgrid arrays (left panel), and sweeps in the y 
direction are only performed in the wide subgrid arrays (right panel). Thus every sweep operation occurs along the largest possible array 
dimension and involves only the arrays kept in the memory of the responsible processor. Communication between processors is only required 
to exchange updated block arrays of the hydrodynamic variables after each sweep is finished. 



